The retrieval of high resolution ground profiles is of great importance for the analysis of geomorphological processes such as flow processes (Mueting, Bookhagen, and Strecker, 2021) and serves as the basis for research on river flow gradient analysis (Scherer et al., 2022) or aboveground biomass estimation (Atmani, Bookhagen, and Smith, 2022). The digital elevation models (DEMs) used for such analyses can be derived through different methods, such as stereophotogrammetry with passively acquired data or active methods like radar (radio detection and ranging) and lidar (light detection and ranging). Lidar derived products in particular can reach high spatial resolutions (Neuenschwander and Pitts, 2019) and are usually acquired from the air or space for DEM generation.

One such spaceborne mission that operates with lidar is ICESat-2 with an identically-named satellite that has the “Advanced Topographic Laser Altimeter System” (ATLAS) instrument onboard. ICESat-2’s objective is to measure and help map the elevation of clouds and ice and land surfaces globally and monitor vegetation height. The laser of ICESat-2 is split into six beams in three pairs, which are approximately 3.3 kilometers apart across-track, the beams of each pair are 90 meters apart. Each pair has a stronger left beam and a weaker right beam with each beam having a footprint of 17 m diameter with a 0.7 m sampling interval (Neuenschwander and Pitts, 2019). ICESat-2’s data products comes in different processing levels with varying features. The level 2 data product “Global Geolocated Photon Data” (ATL03) features all recorded photons, containing information on latitude, longitude, height, surface type and signal confidence. An ICESat-2 product that has global terrain height available is the level 3b “Global Geolocated Photon Data” (ATL08) but it has a fixed downsampled spatial resolution of 100 meters (National Snow and Ice Data Center, 2023). For detailed analyses on landscapes with large changes in relief over small areas, this resolution is too coarse.

Given the necessity for high resolution ground profiles and the promising 70 centimeter sampling interval of ICESat-2 data in the ATL03 product, this raises the following research question:

What are useful filtering steps to extract ground profiles from ICESat-2 ATL03 data and what accuracy can be achieved?

This internship was supervised by Prof. Dr. Bodo Bookhagen.

The study area of this internship consists of the Quebrada del Toro and the Río Capilla catchment west of it, both located in the Toro basin in the Eastern Cordillera of northwestern Argentina in the Central Andes. Figure 1 shows the topography and location of the study area. The Toro basin is located between the Puna-Altiplano Plateau in the west and the foreland in the east, near the city of Salta (Mueting, Bookhagen, and Strecker, 2021; Hilley and Strecker, 2005). The Quebrada del Toro is characterized by high relief with deep incision (Mueting, Bookhagen, and Strecker, 2021). In the northernmost part of the study area, elevation is 3,500 m a.s.l. in the mountains and 2,700 m in the valleys. The central part of the study site shows the highest elevation with up to 4,180 m a.s.l. and the narrow, steep valley of Río Toro that incises the landscape (Mueting, Bookhagen, and Strecker, 2021). The basin outlet in the forelands in the southern part of the study area shows elevations of 1,450 m a.s.l.

ATL03 data can be accessed and downloaded as hdf5 files through the Data Access Tool of the NSIDC (National Snow and Ice Data Center). For this internship, a dataset from 29/05/2022 that goes through the center of the study area was chosen (see Figure 1). The reference ground track of the dataset is 1032, cycle number 15, segment number 8.

Each of the six ICESat-2 beams was extracted and then filtered for signal photons by the confidence flags that ATL03 data provides for every photon: flags 2 to 4 indicate low, medium and high confidence, respectively; photons that were labelled noise or unspecified were not used for analysis. The distribution of photon confidence along the track can be seen in Figure 2.

Validation data was used in the form of DEMs. Three DEMs were used with the addition of a global geoid dataset to correct for different height coordinate reference systems.

The first DEM used was the Quebrada del Toro 3m DEM, which was created in Mueting, Bookhagen, and Strecker, 2021 from tri-stereo SPOT-7 satellite data. The DEM file was derived from Zenodo and has a accuracy of 0.051 ± 1.915 m, but covers a smaller extent than the study area. Its vertical coordinate system is WGS84 (Mueting, Bookhagen, and Strecker, 2021).

Another high-resolution DEM was the Salta Sector 5.1 DEM by the Instituto Geográfico Nacional of Argentina with a spatial resolution of 5 meters. Its vertical reference system is the Sistema de Referencia Vertical Nacional del año 1971 (SRVN71), which is relative to the localized tide gauge measured in Mar del Plata, Argentina (Instituto Geográfico Nacional, 2023).

Additionally the Copernicus DEM in the GLO-30 instance was used. The data covers the global land surface and has a spatial resolution of 30 meters and an absolute vertical accuracy of ground profile photons # assumes residuals are normal distributed for i in range(4): print(f'Loop {i+1}') print(f'Number of photons: {len(xy_loop)}') # calculate linear fit on kN = 150, smooth, remove outliers coeffs, pfit_elev = polyfit_numpy(xy_loop, 150, 1) pfit_elev_smooth = gaussian_filter(pfit_elev, sigma=5) pfit_elev_smooth[np.abs(xy_loop[:,1]-pfit_elev_smooth) > 50] = np.nan print(f'Number of elements in smoothed polyfit: {np.sum(~np.isnan(pfit_elev_smooth))}') # suppress 'All-NaN slice encountered' or 'Mean of empty slice' RuntimeWarnings that will occur with warnings.catch_warnings(): warnings.simplefilter("ignore", category=RuntimeWarning) # statistics on photon residuals residuals = xy_loop[:,1] - pfit_elev_smooth print(f'Max residual at index {np.nanargmax(residuals)} with {np.nanmax(residuals):.2f} m') print(f'Min residual at index {np.nanargmin(residuals)} with {np.nanmin(residuals):.2f} m') print(f'StDev of residuals: {np.nanstd(residuals):.2f}') # split photons and residuals into 30 m bins split_at = xy_loop[:, 0].searchsorted(np.arange(0, int(max(xy_loop[:,0])), 30)) xy_bins = np.split(xy_loop, split_at) residuals_bins = np.split(residuals, split_at) # empty list to catch filtered photon bins new_xy_bins = [] for j in range(len(xy_bins)): # retrieve indices where residuals are in 20-80th residual percentile filtered_ind, = np.where((residuals_bins[j] > np.nanpercentile(residuals_bins[j], 20)) & (residuals_bins[j]